import geopandas as pd
import contextily as ctx # Used for contextual basemaps
import matplotlib.pyplot as plt
from geocube.api.core import make_geocube # Used for rasterizing
import os
import shapely
import imageio # Used for making animated GIFs
import numpy as np
from IPython.display import Image
from osgeo import gdal # Raster operations
import zipfile
import rasterio
import rasterio.merge
import rasterio.plot
import rasterio.warp
plt.rcParams['figure.figsize'] = (20, 20)
os.listdir("input")
/usr/local/lib/python3.8/dist-packages/geopandas/_compat.py:106: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.0-CAPI-1.16.2). Conversions between both will be slow. warnings.warn(
['lds-nz-road-centrelines-topo-150k-FGDB.zip', 'lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip', 'statsnzregional-council-2021-clipped-generalised-FGDB.zip', 'statsnz2018-census-usually-resident-population-and-age-groups-by-te-FGDB.zip', 'Subnational population projections 2018-2048.xls', 'lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip']
First, read regional council bounds. This geometry will be used to clip NZ-wide datasets to just the region of interest, Auckland
%%time
REGC = pd.read_file("input/statsnzregional-council-2021-clipped-generalised-FGDB.zip!regional-council-2021-clipped-generalised.gdb")
AKL = REGC[REGC.REGC2021_V1_00_NAME == "Auckland Region"].copy()
# Filter out islands
AKL["geometry"] = max(AKL.geometry.explode(), key=lambda a: a.area)
# Coordinate reference system (projection)
print(AKL.crs)
# Simplify geometry to speed up clip operations
AKL = AKL.simplify(1000).buffer(1000)
ax = AKL.to_crs(epsg=3857).boundary.plot()
ax.set_title("Auckland Region clip extent")
ctx.add_basemap(ax)
epsg:2193 CPU times: user 1.75 s, sys: 986 ms, total: 2.73 s Wall time: 14.6 s
Load the LRIS Land Cover Database (downloaded in GDB format from https://lris.scinfo.org.nz/layer/104400-lcdb-v50-land-cover-database-version-50-mainland-new-zealand/)
%%time
df = pd.read_file("zip://input/lris-lcdb-v50-land-cover-database-version-50-mainland-new-zealand-FGDB.zip!lcdb-v50-land-cover-database-version-50-mainland-new-zealand.gdb")
CPU times: user 1min 46s, sys: 3.2 s, total: 1min 49s Wall time: 1min 49s
print(df.columns)
print(df.crs)
display(df.sample(5))
Index(['Name_2018', 'Name_2012', 'Name_2008', 'Name_2001', 'Name_1996',
'Class_2018', 'Class_2012', 'Class_2008', 'Class_2001', 'Class_1996',
'Wetland_18', 'Wetland_12', 'Wetland_08', 'Wetland_01', 'Wetland_96',
'Onshore_18', 'Onshore_12', 'Onshore_08', 'Onshore_01', 'Onshore_96',
'EditAuthor', 'EditDate', 'LCDB_UID', 'geometry'],
dtype='object')
epsg:2193
| Name_2018 | Name_2012 | Name_2008 | Name_2001 | Name_1996 | Class_2018 | Class_2012 | Class_2008 | Class_2001 | Class_1996 | ... | Wetland_96 | Onshore_18 | Onshore_12 | Onshore_08 | Onshore_01 | Onshore_96 | EditAuthor | EditDate | LCDB_UID | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 174810 | Indigenous Forest | Indigenous Forest | Indigenous Forest | Indigenous Forest | Indigenous Forest | 69 | 69 | 69 | 69 | 69 | ... | no | yes | yes | yes | yes | yes | Terralink | 2004-06-30T00:00:00 | lcdb2000180191 | MULTIPOLYGON (((1271800.779 4844357.101, 12717... |
| 4423 | Gorse and/or Broom | Gorse and/or Broom | Gorse and/or Broom | Gorse and/or Broom | Gorse and/or Broom | 51 | 51 | 51 | 51 | 51 | ... | no | yes | yes | yes | yes | yes | Terralink | 2004-06-30T00:00:00 | lcdb2000105107 | MULTIPOLYGON (((1528090.924 5194424.609, 15280... |
| 430719 | Lake or Pond | Lake or Pond | Lake or Pond | Lake or Pond | Lake or Pond | 20 | 20 | 20 | 20 | 20 | ... | no | yes | yes | yes | yes | yes | Terralink | 2004-06-30T00:00:00 | lcdb2000045836 | MULTIPOLYGON (((1261101.499 4966013.291, 12611... |
| 461252 | Indigenous Forest | Indigenous Forest | Indigenous Forest | Indigenous Forest | Indigenous Forest | 69 | 69 | 69 | 69 | 69 | ... | no | yes | yes | yes | yes | yes | Terralink | 2004-06-30T00:00:00 | lcdb1000169706 | MULTIPOLYGON (((1710651.930 6020180.126, 17106... |
| 344602 | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | 52 | 52 | 52 | 52 | 52 | ... | yes | yes | yes | yes | yes | yes | Landcare Research | 2014-06-30T00:00:00 | lcdb1000098963 | MULTIPOLYGON (((1612278.405 6152038.608, 16122... |
5 rows × 24 columns
%%time
df = pd.clip(df, AKL)
CPU times: user 46.4 s, sys: 8.55 ms, total: 46.4 s Wall time: 46.4 s
df.sample(5)
| Name_2018 | Name_2012 | Name_2008 | Name_2001 | Name_1996 | Class_2018 | Class_2012 | Class_2008 | Class_2001 | Class_1996 | ... | Wetland_96 | Onshore_18 | Onshore_12 | Onshore_08 | Onshore_01 | Onshore_96 | EditAuthor | EditDate | LCDB_UID | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 173780 | Urban Parkland/Open Space | Urban Parkland/Open Space | Urban Parkland/Open Space | Urban Parkland/Open Space | Urban Parkland/Open Space | 2 | 2 | 2 | 2 | 2 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2004-06-30T00:00:00 | lcdb1000006473 | POLYGON ((1757451.024 5944934.676, 1757463.730... |
| 57529 | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | 52 | 52 | 52 | 52 | 52 | ... | no | yes | yes | yes | yes | yes | Terralink | 2004-06-30T00:00:00 | lcdb1000093800 | POLYGON ((1733473.489 5921727.296, 1733459.367... |
| 408632 | Low Producing Grassland | Low Producing Grassland | Exotic Forest | Exotic Forest | Exotic Forest | 41 | 41 | 71 | 71 | 71 | ... | no | yes | yes | yes | yes | yes | Landcare Research | 2019-12-01T00:00:00 | lcdb1000415656 | POLYGON ((1778637.742 5904919.429, 1778637.742... |
| 394053 | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | Manuka and/or Kanuka | 52 | 52 | 52 | 52 | 52 | ... | no | yes | yes | yes | yes | yes | Terralink | 2004-06-30T00:00:00 | lcdb1000094687 | POLYGON ((1735422.343 5948607.284, 1735382.309... |
| 507426 | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | Herbaceous Saline Vegetation | 46 | 46 | 46 | 46 | 46 | ... | yes | yes | yes | yes | yes | yes | Terralink | 2004-06-30T00:00:00 | lcdb1000066435 | POLYGON ((1749641.378 5917711.857, 1749687.537... |
5 rows × 24 columns
df.Name_2018.value_counts()
Exotic Forest 3981 Indigenous Forest 3673 Manuka and/or Kanuka 2282 Broadleaved Indigenous Hardwoods 1788 Built-up Area (settlement) 1350 High Producing Exotic Grassland 1326 Mangrove 1151 Urban Parkland/Open Space 1099 Estuarine Open Water 441 Orchard, Vineyard or Other Perennial Crop 436 Short-rotation Cropland 362 Lake or Pond 326 Herbaceous Saline Vegetation 303 Low Producing Grassland 291 Gorse and/or Broom 287 Forest - Harvested 266 Sand or Gravel 252 Deciduous Hardwoods 201 Surface Mine or Dump 132 Mixed Exotic Shrubland 120 Herbaceous Freshwater Vegetation 118 Transport Infrastructure 107 River 15 Flaxland 9 Gravel or Rock 9 Fernland 2 Matagouri or Grey Scrub 1 Name: Name_2018, dtype: int64
These classes are far too detailed - simplify to just Urban, Vegetation, Water, Other
def simplify_classes(code):
if code in [1, 2, 5]:
return 1, "Urban"
elif code in [68,69,71]:
return 2, "Vegetation"
elif code in [0,20,21,22,45,46]:
return 3, "Water"
else:
return 4, "Other"
summary = []
years = [1996, 2001, 2008, 2012, 2018]
for year in years:
print(year)
class_year = f"Class_{year}"
df[class_year + "_simplified_code"] = df[class_year].apply(lambda c: simplify_classes(c)[0])
df[class_year + "_simplified_name"] = df[class_year].apply(lambda c: simplify_classes(c)[1])
summary.append(df[class_year + "_simplified_name"].value_counts())
1996 2001 2008 2012 2018
pd.GeoDataFrame(summary).plot.area()
<AxesSubplot:>
%%capture
# %%capture suppresses output
if not os.path.isfile("land_use.gif"):
ims = []
years = [1996, 2001, 2008, 2012, 2018]
for year in years:
ax = df.plot(column=f'Class_{year}_simplified_name', legend=True)
ax.set_title(year)
ax.figure.tight_layout()
canvas = ax.figure.canvas
canvas.draw() # draw the canvas, cache the renderer
image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
image = image.reshape(canvas.get_width_height()[::-1] + (3,))
ims.append(image)
imageio.mimsave("land_use.gif", ims, fps=1)
with open('land_use.gif','rb') as file:
display(Image(file.read()))
cols = [f"Class_{year}_simplified_code" for year in years]
cols
['Class_1996_simplified_code', 'Class_2001_simplified_code', 'Class_2008_simplified_code', 'Class_2012_simplified_code', 'Class_2018_simplified_code']
%%time
geocube = make_geocube(
vector_data=df,
output_crs="epsg:2193",
measurements=cols,
resolution=(-100, 100),
fill=0, # NaNs, like offshore areas, will be 0
)
geocube
CPU times: user 18.4 s, sys: 9.83 ms, total: 18.4 s Wall time: 18.5 s
<xarray.Dataset>
Dimensions: (x: 1001, y: 1320)
Coordinates:
* y (y) float64 6.002e+06 6.002e+06 ... 5.87e+06
* x (x) float64 1.704e+06 1.704e+06 ... 1.804e+06
spatial_ref int64 0
Data variables:
Class_1996_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Class_2001_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Class_2008_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Class_2012_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Class_2018_simplified_code (y, x) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0
Attributes:
grid_mapping: spatial_refarray([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])
array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])
array(0)
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])geocube.Class_2018_simplified_code.plot()
<matplotlib.collections.QuadMesh at 0x7f6c907ea910>
for year in years:
print(year)
outfile = f"output/land_use_{year}.tif"
if not os.path.isfile(outfile):
geocube[f"Class_{year}_simplified_code"].rio.to_raster(outfile, dtype=np.byte) # Use np.byte for smaller output filesize
1996 2001 2008 2012 2018
%%time
pop = pd.read_file("input/statsnz2018-census-usually-resident-population-and-age-groups-by-te-FGDB.zip!2018-census-usually-resident-population-and-age-groups-by-te.gdb")
pop = pd.clip(pop, AKL)
pop.head()
CPU times: user 4.16 s, sys: 0 ns, total: 4.16 s Wall time: 4.15 s
| TALB2018_V1_00 | TALB2018_V1_00_NAME | LAND_AREA_SQ_KM | AREA_SQ_KM | Age_life_cycle_groups_15_29_years_2006 | Age_life_cycle_groups_30_64_years_2006 | Age_life_cycle_groups_65_years_and_over_2006 | Age_life_cycle_groups_total_2006 | Age_life_cycle_groups_Under_15_years_2006 | Age_life_cycle_groups_15_29_years_2013 | ... | Pop_Total_2006 | Pop_Total_2013 | Pop_Total_2018 | Age_life_cycle_groups_15_29_years_2013_2018_percentChange | Age_life_cycle_groups_30_64_years_2013_2018_percentChange | Age_life_cycle_groups_65_years_and_over_2013_2018_percentChan | Age_life_cycle_groups_total_2013_2018_percentChange | Age_life_cycle_groups_Under_15_years_2013_2018_percentChange | Shape_Length | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 2 | 00300 | Kaipara District | 3108.710953 | 3114.579850 | 2631 | 8568 | 2715 | 18135 | 4215 | 2670 | ... | 18135 | 18963 | 22869 | 25.6 | 15.2 | 33.8 | 20.6 | 16.4 | 902974.494161 | MULTIPOLYGON (((1718576.593 5980699.195, 17185... |
| 4 | 01200 | Hauraki District | 1270.087068 | 1270.087068 | 2568 | 8274 | 3051 | 17856 | 3963 | 2634 | ... | 17856 | 17808 | 20022 | 12.3 | 10.2 | 21.3 | 12.4 | 7.5 | 311223.430113 | POLYGON ((1801726.041 5899678.947, 1803788.474... |
| 5 | 01300 | Waikato District | 4403.282847 | 4516.062518 | 9888 | 27450 | 5580 | 57588 | 14670 | 10695 | ... | 57588 | 63378 | 75618 | 23.8 | 17.2 | 26.9 | 19.3 | 16.6 | 491881.373163 | POLYGON ((1747419.052 5870470.171, 1747413.783... |
| 29 | 07606 | Henderson-Massey Local Board Area | 53.219988 | 64.608373 | 21441 | 43914 | 9264 | 98790 | 24174 | 23274 | ... | 98790 | 107685 | 118422 | 12.7 | 10.6 | 8.3 | 10.0 | 7.0 | 47587.459227 | MULTIPOLYGON (((1748774.713 5920936.996, 17483... |
| 30 | 07607 | Waitakere Ranges Local Board Area | 305.732502 | 943.004892 | 8343 | 22752 | 3372 | 45498 | 11031 | 8763 | ... | 45498 | 48399 | 52095 | 13.4 | 5.6 | 18.6 | 7.6 | 3.0 | 148724.962888 | POLYGON ((1742888.086 5900837.328, 1742806.519... |
5 rows × 29 columns
pop.columns
Index(['TALB2018_V1_00', 'TALB2018_V1_00_NAME', 'LAND_AREA_SQ_KM',
'AREA_SQ_KM', 'Age_life_cycle_groups_15_29_years_2006',
'Age_life_cycle_groups_30_64_years_2006',
'Age_life_cycle_groups_65_years_and_over_2006',
'Age_life_cycle_groups_total_2006',
'Age_life_cycle_groups_Under_15_years_2006',
'Age_life_cycle_groups_15_29_years_2013',
'Age_life_cycle_groups_30_64_years_2013',
'Age_life_cycle_groups_65_years_and_over_2013',
'Age_life_cycle_groups_total_2013',
'Age_life_cycle_groups_Under_15_years_2013',
'Age_life_cycle_groups_15_29_years_2018',
'Age_life_cycle_groups_30_64_years_2018',
'Age_life_cycle_groups_65_years_and_over_2018',
'Age_life_cycle_groups_total_2018',
'Age_life_cycle_groups_Under_15_years_2018', 'Pop_Total_2006',
'Pop_Total_2013', 'Pop_Total_2018',
'Age_life_cycle_groups_15_29_years_2013_2018_percentChange',
'Age_life_cycle_groups_30_64_years_2013_2018_percentChange',
'Age_life_cycle_groups_65_years_and_over_2013_2018_percentChan',
'Age_life_cycle_groups_total_2013_2018_percentChange',
'Age_life_cycle_groups_Under_15_years_2013_2018_percentChange',
'Shape_Length', 'geometry'],
dtype='object')
cols = ["Pop_Total_2006", "Pop_Total_2013", "Pop_Total_2018"]
d = pop[cols].to_numpy()
d.min(), d.mean(), d.max()
(39, 63916.25, 140970)
%%capture
# %%capture suppresses output
if not os.path.isfile("pop.gif"):
ims = []
for col in cols:
ax = pop.plot(column=col, legend=True, vmin=0, vmax=d.max(), cmap='OrRd')
ax.set_title(col)
ax.figure.tight_layout()
canvas = ax.figure.canvas
canvas.draw() # draw the canvas, cache the renderer
image = np.frombuffer(canvas.tostring_rgb(), dtype='uint8')
image = image.reshape(canvas.get_width_height()[::-1] + (3,))
ims.append(image)
imageio.mimsave("pop.gif", ims, fps=1)
with open('pop.gif','rb') as file:
display(Image(file.read()))
%%time
popcube = make_geocube(
vector_data=pop,
measurements=["Pop_Total_2006", "Pop_Total_2013", "Pop_Total_2018"],
like=geocube, # Ensures dimensions match
fill=0 # NaNs, like offshore areas, will be 0
)
popcube
CPU times: user 271 ms, sys: 19.5 ms, total: 290 ms Wall time: 289 ms
<xarray.Dataset>
Dimensions: (x: 1001, y: 1320)
Coordinates:
* y (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
* x (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
spatial_ref int64 0
Data variables:
Pop_Total_2006 (y, x) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Pop_Total_2013 (y, x) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Pop_Total_2018 (y, x) float64 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
grid_mapping: spatial_refarray([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])
array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])
array(0)
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])for col in cols:
outfile = f"output/{col}.tif"
if not os.path.isfile(outfile):
# byte max value is 255, and we have larger values than that here. uint32 max value is 4294967295
popcube[col].rio.to_raster(outfile, dtype=np.uint32)
%%time
roads = pd.read_file("input/lds-nz-road-centrelines-topo-150k-FGDB.zip!nz-road-centrelines-topo-150k.gdb")
CPU times: user 8.6 s, sys: 20 ms, total: 8.62 s Wall time: 8.62 s
%%time
akl_roads = pd.clip(roads, AKL)
CPU times: user 25.2 s, sys: 0 ns, total: 25.2 s Wall time: 25.2 s
# If a road has a highway number (hway_num not None), it's a highway/motorway
mway = akl_roads[~akl_roads.hway_num.isna()].copy()
mway
| t50_fid | name_ascii | macronated | name | hway_num | rna_sufi | lane_count | way_count | status | surface | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 512 | 100120610 | KAIPARA COAST HIGHWAY | N | KAIPARA COAST HIGHWAY | 16 | 3007739 | 2 | None | None | sealed | LINESTRING (1732000.000 5944172.070, 1732048.5... |
| 2933 | 3198057 | STATE HIGHWAY 1 | N | STATE HIGHWAY 1 | 1 | 3027695 | 2 | None | None | sealed | LINESTRING (1748581.508 5968975.145, 1748558.4... |
| 2934 | 3198059 | STATE HIGHWAY 1 | N | STATE HIGHWAY 1 | 1 | 3027695 | 2 | None | None | sealed | LINESTRING (1748171.047 5971284.152, 1748129.9... |
| 3320 | 3200754 | PAERATA ROAD | N | PAERATA ROAD | 22 | 3000260 | 2 | None | None | sealed | LINESTRING (1767236.112 5888088.508, 1767244.3... |
| 3324 | 3200792 | UPPER HARBOUR MOTORWAY | N | UPPER HARBOUR MOTORWAY | 18 | 3047073 | 4 | None | None | sealed | LINESTRING (1747954.314 5927269.837, 1747970.0... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 138240 | 100048291 | AUCKLAND-WAIWERA MOTORWAY | N | AUCKLAND-WAIWERA MOTORWAY | 1 | 3067966 | 7 | None | None | sealed | LINESTRING (1755881.018 5922863.734, 1755886.4... |
| 138301 | 100048432 | AUCKLAND-HAMILTON MOTORWAY | N | AUCKLAND-HAMILTON MOTORWAY | 1 | 3017109 | 1 | None | None | sealed | LINESTRING (1765115.647 5909916.697, 1765092.7... |
| 138337 | 100048532 | STATE HIGHWAY 1 | N | STATE HIGHWAY 1 | 1 | 3027695 | 4 | None | None | sealed | LINESTRING (1748892.089 5949596.727, 1748892.0... |
| 138369 | 100048589 | PORT ALBERT ROAD | N | PORT ALBERT ROAD | 16 | 3013274 | 2 | None | None | sealed | LINESTRING (1734173.019 5980575.187, 1734175.6... |
| 138680 | 100118365 | SOUTH-WESTERN MOTORWAY | N | SOUTH-WESTERN MOTORWAY | 20 | 3018532 | 4 | None | None | sealed | LINESTRING (1760066.252 5908184.133, 1760043.9... |
426 rows × 11 columns
mway.name.value_counts().head(50).plot(kind="barh").invert_yaxis()
mway.hway_num.value_counts().head(50).plot(kind="barh").invert_yaxis()
ax = mway.to_crs(epsg=3857).plot()
ax.set_title("Auckland Region motorways")
ctx.add_basemap(ax)
%%time
mway_cube = make_geocube(
vector_data=mway,
measurements=["lane_count"],
like=geocube, # Ensures dimensions match
fill=0, # 0 works fine here, as every mway has at least one lane
)
mway_cube
CPU times: user 259 ms, sys: 0 ns, total: 259 ms Wall time: 256 ms
<xarray.Dataset>
Dimensions: (x: 1001, y: 1320)
Coordinates:
* y (y) float64 6.002e+06 6.002e+06 ... 5.871e+06 5.87e+06
* x (x) float64 1.704e+06 1.704e+06 ... 1.804e+06 1.804e+06
spatial_ref int64 0
Data variables:
lane_count (y, x) float64 0.0 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes:
grid_mapping: spatial_refarray([6002350., 6002250., 6002150., ..., 5870650., 5870550., 5870450.])
array([1703950., 1704050., 1704150., ..., 1803750., 1803850., 1803950.])
array(0)
array([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]])mway_cube.lane_count.plot()
outfile = "output/mway.tif"
if not os.path.isfile(outfile):
mway_cube.lane_count.rio.to_raster(outfile, dtype=np.byte)
src_ds = gdal.Open("output/mway.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/mway_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
src_ds.RasterXSize, src_ds.RasterYSize, 1,
gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
mway_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(mway_dist.shape)
plt.imshow(mway_dist)
plt.title("Distance from motorways in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")
(1320, 1001)
Text(0.5, 1.0, 'Distance (m)')
%%time
cbd = pop2018[pop2018.MB2020_V2_00 == "0433501"].copy()
cbd.geometry = cbd.geometry.buffer(1000)
cbd_cube = make_geocube(
vector_data=cbd,
like=geocube, # Ensures dimensions match
fill=0
)
outfile = "output/cbd.tif"
if not os.path.isfile(outfile):
cbd_cube.General_Electoral_Population.rio.to_raster(outfile, dtype=np.byte)
src_ds = gdal.Open("output/cbd.tif")
srcband = src_ds.GetRasterBand(1)
dst_filename = "output/cbd_dist.tif"
drv = gdal.GetDriverByName('GTiff')
dst_ds = drv.Create( dst_filename,
src_ds.RasterXSize, src_ds.RasterYSize, 1,
gdal.GetDataTypeByName('UInt16'))
dst_ds.SetGeoTransform( src_ds.GetGeoTransform() )
dst_ds.SetProjection( src_ds.GetProjectionRef() )
dstband = dst_ds.GetRasterBand(1)
prox = gdal.ComputeProximity(srcband,dstband,["DISTUNITS=GEO"]) # Encoded value is distance from motorway in meters
# Garbage collection of this variable flushes write
dst_ds = None
dst_ds = gdal.Open(dst_filename)
cbd_dist = np.array(dst_ds.GetRasterBand(1).ReadAsArray())
print(cbd_dist.shape)
plt.imshow(cbd_dist)
plt.title("Distance from CBD in Auckland")
cb = plt.colorbar()
cb.ax.set_title("Distance (m)")
(1320, 1001) CPU times: user 441 ms, sys: 10 ms, total: 451 ms Wall time: 446 ms
Text(0.5, 1.0, 'Distance (m)')
bounds = AKL.total_bounds.tolist()
bounds
[1703081.9789640256, 5870396.320936217, 1804839.668875325, 6002367.198185163]
zf = zipfile.ZipFile('input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip')
tiles = [file for file in zf.namelist() if file.endswith(".tif")]
tiles
['EJ.tif', 'DM.tif', 'EL.tif', 'DL.tif', 'DJ.tif', 'FK.tif', 'DK.tif', 'EK.tif', 'FL.tif']
tile_datasets = [rasterio.open(f'zip://input/lds-nz-8m-digital-elevation-model-2012-GTiff-auckland-region.zip!{tile}') for tile in tiles]
DEM, transform = rasterio.merge.merge(tile_datasets, bounds = bounds, res = (100,100), dtype=np.int16)
print(np.nanmin(DEM), np.nanmean(DEM), np.nanmax(DEM), DEM.shape)
rasterio.plot.show(np.where(DEM>=0, DEM, np.nan), cmap='terrain')
-32767 -18316.517014943143 697 (1, 1320, 1018)
<AxesSubplot:>
cbd_dist = rasterio.open("output/cbd_dist.tif")
cbd_dist.read().shape
(1, 1320, 1001)
NODATA = -32767
warped_DEM = np.full_like(cbd_dist.read(), NODATA, dtype=np.int16)
rasterio.warp.reproject(DEM,
warped_DEM,
src_transform=transform,
src_nodata=NODATA,
dst_nodata=NODATA,
src_crs=tile_datasets[0].crs,
dst_crs=cbd_dist.crs,
dst_transform=cbd_dist.transform
)
print(DEM.shape, warped_DEM.shape)
(1, 1320, 1018) (1, 1320, 1001)
meta = tile_datasets[0].meta
print(meta)
meta.update({
"dtype": "int16",
"height": cbd_dist.height,
"width": cbd_dist.width,
"transform": cbd_dist.transform
})
print(meta)
outfile = "output/dem.tif"
if not os.path.isfile(outfile):
with rasterio.open(outfile, "w", **meta) as dest:
dest.write(warped_DEM)
{'driver': 'GTiff', 'dtype': 'float32', 'nodata': -32767.0, 'width': 3028, 'height': 8192, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(8.0, 0.0, 1679712.0,
0.0, -8.0, 5963776.0)}
{'driver': 'GTiff', 'dtype': 'int16', 'nodata': -32767.0, 'width': 1001, 'height': 1320, 'count': 1, 'crs': CRS.from_epsg(2193), 'transform': Affine(100.0, 0.0, 1703900.0,
0.0, -100.0, 6002400.0)}
gdal.DEMProcessing('output/slope.tif', outfile, 'slope')
slope = rasterio.open("output/slope.tif").read()
print(np.nanmin(slope), np.nanmean(slope), np.nanmax(slope), slope.shape)
rasterio.plot.show(np.where(slope>=0, slope, np.nan), cmap='terrain')
-9999.0 -5778.3354 52.492313 (1, 1320, 1001)
<AxesSubplot:>
rasters = [rasterio.open(f"output/{f}") for f in os.listdir("output") if f.endswith(".tif")]
params = set([(r.shape, r.res, r.crs, r.count, r.bounds) for r in rasters])
print(params)
# Assert all rasters have the same shape, pixel size, CRS, number of bands, and bounds
assert len(params) == 1
{((1320, 1001), (100.0, 100.0), CRS.from_epsg(2193), 1, BoundingBox(left=1703900.0, bottom=5870400.0, right=1804000.0, top=6002400.0))}
!ls -Ggh output
total 27M -rw-r--r-- 1 1.3M Apr 16 14:44 cbd.tif -rw-r--r-- 1 2.6M Apr 16 16:09 cbd_dist.tif -rw-r--r-- 1 2.6M Apr 16 16:46 dem.tif -rw-r--r-- 1 1.3M Apr 16 09:50 land_use_1996.tif -rw-r--r-- 1 1.3M Apr 16 09:50 land_use_2001.tif -rw-r--r-- 1 1.3M Apr 16 09:50 land_use_2008.tif -rw-r--r-- 1 1.3M Apr 16 09:50 land_use_2012.tif -rw-r--r-- 1 1.3M Apr 16 09:50 land_use_2018.tif -rw-r--r-- 1 1.3M Apr 16 09:50 mway.tif -rw-r--r-- 1 2.6M Apr 16 16:09 mway_dist.tif -rw-r--r-- 1 2.6M Apr 16 09:50 pop2013.tif -rw-r--r-- 1 2.6M Apr 16 09:50 pop2018.tif -rw-r--r-- 1 5.1M Apr 16 16:52 slope.tif